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ANALYSIS OF COMPOSITE SKIN-STIFFENER DEBOND SPECIMENS USING 
A SHELL/3D MODELING TECHNIQUE AND SUBMODELING 

Ronald Krueger 1 and Pierre J. Minguet 2 


ABSTRACT 

The application of a shell/3D modeling technique for the simulation of skin/stringer 
debond in a specimen subjected to tension and three-point bending was studied. The 
global structure was modeled with shell elements. A local three-dimensional model, 
extending to about three specimen thicknesses on either side of the delamination front 
was used to model the details of the damaged section. Computed total strain energy 
release rates and mixed-mode ratios obtained from shell/3 D simulations were in good 
agreement with results obtained from full solid models. The good correlation of the 
results demonstrated the effectiveness of the shell/3D modeling technique for the 
investigation of skin/stiffener separation due to delamination in the adherents. In 
addition, the application of the submodeling technique for the simulation of skin/stringer 
debond was also studied. Global models made of shell elements and solid elements were 
studied. Solid elements were used for local submodels, which extended between three 
and six specimen thicknesses on either side of the delamination front to model the details 
of the damaged section. Computed total strain energy release rates and mixed-mode 
ratios obtained from the simulations using the submodeling technique were not in 
agreement with results obtained from full solid models. 


1. BACKGROUND 

Many composite components in aerospace structures are made of flat or curved panels 
with co-cured or adhesively bonded frames and stiffeners. Recent studies focused on the 
investigation of the debonding mechanism and included testing of skin/stiffener panels and 
failure analysis using shell models [1-4]. Over the last decade a consistent step-wise approach 
has been developed which uses experiments to detect the failure mechanism, computational 
stress analysis to determine the location of first matrix cracking and computational fracture 
mechanics to investigate the potential for delamination growth. Testing of skin gage stiffened 
panels designed for pressurized aircraft fuselage has shown that bond failure at the tip of the 
frame flange is an important and very likely failure mode [5]. Comparatively simple specimens 
consisting of a stringer flange bonded onto a skin have been developed to study skin/stiffener 
debonding [6-8]. The failure that initiates at the tip of the flange in these specimens is nearly 
identical to the failure observed in the full-scale panels and the frame pull-off specimens [7, 9, 
10 ]. 

During previous analyses of skin/stringer debond specimens two-dimensional finite 
element models as shown in Figure 1(a) were used because modeling and computational times 
are shorter compared to full three-dimensional analysis, especially if many different 
configurations have to be analyzed during the initial design phase of a structure [7, 9, 10]. 
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Results from a study where energy release rates from two-dimensional analyses were compared 
to data from full three-dimensional simulations indicated that plane stress and plane strain 
models yield upper and lower-bound results [11]. Accurate predictions, however, require the use 
of three-dimensional models as shown in Figure 1(b). For detailed modeling and analysis of the 
damage observed during experiments, the shell/3D modeling technique - as shown in Figure 1(c) 
- is likely to reduce the modeling effort and save computer time because only a relatively small 
section of interest needs to modeled with solid elements. The technique combines the accuracy 
of the full three-dimensional solution with the computational efficiency of a plate or shell finite 
element model and has been demonstrated for various applications [12-16]. In future applications 
the shell/3D modeling technique therefore is expected to be used for structural components such 
as stringer stiffened panels as outlined in Figure 1(d), where the savings in modeling effort and 
computer time will be even more significant [17]. 

The first objective of this study is to demonstrate the use of the shell/3 D modeling 
technique for the simulation of skin/stringer debond specimens subjected to tension and three- 
point bending. A local three-dimensional model, extending to about three specimen thicknesses 
on either side of the delamination front was used. The transition from shell to solid model was 
accomplished using a new surface based method in ABAQUS® 6.3 to connect the meshed 
regions [18]. Delaminations of various lengths were discretely modeled at the locations where 
delaminations were observed during the experiments. This study complements an investigation 
of the skin/stringer debond specimen subjected to bending where the connection was based on 
individual multi point constraints at the shell/solid interface [16, 19]. 

The second objective of this study is to examine the use of the submodeling technique in 
ABAQUS® for the simulation of skin/stringer debond specimen subjected to tension and three- 
point bending. Global models made of shell elements and solid elements were studied. Local 
submodels were made of solid elements and extended between three and six specimen 
thicknesses on either side of the delamination front. 

For both modeling techniques, mixed mode strain energy release rate distributions were 
computed across the width of the specimens using the virtual crack closure technique (VCCT) 
[20, 21], The results were compared to mixed mode strain energy release rates obtained from 
computations where the entire specimen had been modeled with solid elements [11, 16, 22], 


2. MATERIALS AND TESTING 

The current study uses experimental measurements of skin-stiffener debond specimens as 
reference. The specimens consisted of a bonded skin and flange assembly as illustrated in 
Figure 2(a) and described in reference [10]. An IM7/8552 graphite/epoxy system was used for 
both the skin and flange. The skin was made of prepreg tape with a measured average ply 
thickness of h =0.148 mm and had a [45/-45/0/-45/45/90/90/-45/45/0/45/-45] lay-up. The flange 
was made of a plain-weave fabric with a thickness of h =0.212 mm. The flange lay-up was 
[45/0/45/0/45/0/45/0/45] f, where the subscript “f” denotes fabric, “0” represents a 0°-90° fabric 
ply and “45” represents a 0°-90° fabric ply rotated by 45° [23]. The measured bondline thickness 
averaged 0.178mm. Specimens were 25.4-mm wide and 177.8-mm long. The properties of the 
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graphite/epoxy material and the adhesive were measured at Boeing and are summarized in 
Table 1 

Four quasi-static tension tests were performed in a servohydraulic load frame in 
displacement control at a stroke rate of 0.4 mm/min. The specimens were mounted in hydraulic 
grips with a gage length of 101.6 mm as shown in Figure 2(b). The value of the damage-onset 
load was averaged from four tests and determined to be P =17.8 kN with a coefficient of 
variation of 8.9%. Three-point bending tests were performed in a servohydraulic load frame at a 
monotonic rate of 1.52 mm/min. The bottom support had a span of 101.6 mm as shown in 
Figure 2(c). The value of the applied load at damage-onset was averaged from four tests and 
determined to be Q =427.6 N with a coefficient of variation of 12.8%. The tests were terminated 
when the flange debonded from the skin. Damage was documented from photographs of the 
polished specimen edges at each of the four flange comers identified in Figure 2(a). Comers 1 
and 4 and comers 2 and 3 showed similar damage patterns for both tests. At corners 1 and 4, a 
delamination running in the bondline (delamination Al) initiated from a matrix crack in the 
adhesive pocket as shown in Figure 3(a). In two of the quasi-static cases, a second delamination 
(A2) was observed below the first, in the top — 45°/0° skin ply interface. The damage at comers 2 
and 3, however, formed first and consisted of a matrix crack in the 45° skin surface ply and a 
delamination at the +457-45° interface as shown in Figure 3(b). Therefore, this damage pattern 
has been the focus of the current and previous analyses [10, 11, 16, 21, 22], 


3. FINITE ELEMENT MODELING 

The complex nature of the failure observed during the experiments, where the 
delamination changed across the specimen width from a delamination running at the skin surface 
457-45° layer interface to a delamination propagating in the bondline (see Figure 3), suggests 
the need for a three-dimensional model. The current study presents an intermediate step where 
the local three-dimensional models were created by extruding two-dimensional models across 
the width. The fact that the delamination changed across the specimen width from a delamination 
running at the skin surface 457-45° layer interface to a delamination propagating in the bondline 
above, however, is still not accounted for in this model. Nevertheless, the three-dimensional 
model takes width effects into account and therefore provides additional insight into the 
significance of edge effects. 


3.1. Shell/3D Finite Element Model 

The first objective of this study is to demonstrate a shell/3D modeling technique for the 
investigation of delamination onset from an initial crack in the skin/stringer specimen. The 
shell/3D modeling technique used a local three-dimensional solid finite element model only in 
the immediate vicinity of the delamination front, as shown in the example of a modeled Double 
Cantilever Beam (DCB) specimen in Figure 4(a) [15]. The approach combined the accuracy of 
the full three-dimensional solution with the computational efficiency of a plate or shell finite 
element model. For a local three-dimensional model, extending to a minimum of about three 
specimen thicknesses on either side of the delamination front (d/2h=6.6 7), the computed mixed 
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mode strain energy release rates were within 1% of the results obtained from computations 
where the entire specimen had been modeled with solid elements as shown in Figure 4(b) [15]. 

The shell/3D model of the skin/stringer specimen is shown in Figure 5. The outline in 
Figure 5(a) shows the individual sections and the interface region. The global section was 
modeled with ABAQUS® reduced integrated eight-noded quadrilateral shell elements S8R and 
the local three-dimensional section was modeled with solid C3D20R elements as shown in 
Figure 5(b). The local three-dimensional model extended to about three specimen thicknesses on 
either side of the delamination front as shown in Figure 5(c) since this configuration yielded 
results which were within 1% of the results obtained from full three-dimensional analysis as 
described above. 

The deformed shell/3D model of the specimen subjected to tension loading with load and 
boundary conditions is illustrated in Figure 6(a). The deformation is enlarged to visualize the 
bending of the specimen caused by the load eccentricity in the flange region and the asymmetric 
layup with respect to the neutral axis. The deformed model of the specimen subjected to three- 
point bending with load and boundary conditions is shown in Figure 6(b). Both models have 
identical meshes and differ only in load and boundary conditions. A detail of the modeled 
delaminated region is shown in Figure 6(c). A refined three-dimensional mesh was used in the 
critical area of the 45° skin ply where cracking was observed during the tests as shown in 
Figure 6(c). Outside the refined area, the mesh transitioned rapidly to prevent the three- 
dimensional model from becoming excessively large. The skin plies were grouped into four 
layered elements with 45/-45/0, -45/45/90, 90/-45/45 and 0/45/-45 respectively, thus taking 
advantage of the composite solid element option in ABAQUS® [18]. The fabric layers and the 
resin layer were grouped into five layered elements. In the transition regions several plies were 
modeled by one element with material properties smeared using the rule of mixtures [24], This 
procedure did not calculate the full A-B-D stiffness matrix contributions of the plies [24]; 
however, it appeared suitable for the small transition regions to enforce a reasonable model size. 
The transition from the global shell element model to the local three-dimensional model in the 
vicinity of the delamination front was accomplished by using a new shell to solid coupling 
option offered by ABAQUS® 6.3 to enforce appropriate translations and rotations at the shell- 
solid interface [18]. Using the shell/3D approach the total number of freedoms could be reduced 
by more than 12% compared to the full three-dimensional model of the skin/stringer specimen. 
The efficiency is discussed later. 

The shell/3D model shown in Figure 6, however, was not fine enough in the vicinity of 
the free edges (z=0.0 mm and z=25.4 mm) to accurately model the effect of the free edges on the 
distribution of the energy release rates across the width. Therefore, a shell/3D model as shown in 
Figure 7 was created with twenty elements across the width of the specimen corresponding to the 
full three-dimensional model discussed in reference [16]. The new shell/3D model was improved 
by a refined zone (0.74 mm, five elements) near the free edges (z=0.0 mm and z=25.4 mm) as 
shown in Figures 7(b) and (c). 

The virtual crack closure technique (VCCT) was used to calculate the mode I, II and III 
components of the strain energy release rate for the modeled delamination [20, 21], For the entire 
investigation, the ABAQUS® geometric nonlinear analysis procedure was used. This was done in 
accordance with studies from which reference solutions were taken [11, 16, 22], 
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3.2 Global Finite Element Model and Submodel 

The second objective of this study was to demonstrate the submodeling technique for the 
investigation of delamination onset from an initial crack in the skin/stringer specimen. The 
submodeling technique is generally used to study a local part of a model (submodel) with a 
refined mesh based on interpolation of the solution from an initial, relatively coarse, global 
mesh. The technique is most useful when it is necessary to obtain an accurate, detailed solution 
in a local region and the detailed modeling of that local region has negligible effect on the 
overall solution. One advantage is that different element types can be used in the submodel than 
those used to model the corresponding region in the global model. This allows the use of shell 
elements or solid elements for the global model of the specimen and three-dimensional solid 
elements for the submodel of the delaminated section [18]. 

An overview of the submodeling technique is given in Figure 8 for the case where shell 
elements were used to create the global model. The outline in Figure 8(a) shows the physical 
location of the submodel with respect to the global model of the specimen. To create the global 
model, the entire skin/stinger specimen was modeled with reduced integrated ABAQUS® four- 
noded quadrilateral shell elements S4R or reduced integrated eight-noded quadrilateral shell 
elements S8R as shown in Figure 8(b). The nodes on the submodel boundary where values are 
defined by interpolating the solution from the global model are called “driven nodes”. The 
outline in Figure 8(c) shows the location of the driven nodes on both ends of the submodel. 

The models used for the case where shell elements were used to create the global model 
are shown in Figure 9. The deformed global models of the skin/stringer specimen subjected to 
three-point bending and tension are shown in Figures 9(a) and (b) respectively. The local three- 
dimensional submodel shown in Figure 9(c) was made of solid C3D20R elements. The local 
three-dimensional model extended to about three specimen thicknesses on either side of the 
delamination front and is identical to the local model used in the shell/3D analysis discussed 
above (see Figure 5(c)). This configuration yielded good results during the initial feasibility 
study of the shell/3D technique shown in Figure 4. Therefore, the submodel extending about 
three specimen thicknesses in front and behind the damaged section was used first. In order to 
determine what portion of the structural component needs to be included in a submodel, another 
three submodels as shown in Figure 10 where included in this study. The submodels extended 
about four, five and six specimen thicknesses in front and behind the damaged section. 

For the case where solid elements were used to create the global model, an overview of 
the submodeling technique is given in Figure 11. The outline in Figure 1 1(a) shows the physical 
location of the submodel with respect to the global model of the specimen. To create the global 
model, the entire skin/stinger specimen with was modeled with ABAQUS® reduced integrated 
twenty-noded solid elements C3D20R as shown in Figure 11(b). The outline in Figure 11(c) 
shows the location of the nodes driven by the solution from the global model on both ends of the 
submodel. The submodels are identical to those used in the analysis discussed above. 
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4. FINITE ELEMENT ANALYSIS OF THE SKIN/STRINGER SPECIMEN 

4.1. Analysis of the Skin/Stringer Specimen Subjected to Tension Load Using Shell/3D 

Modeling Technique 

4.1.1 Global Response of the Specimen Subjected to Tension Loading 

Initially, the global response of the specimens was computed. The load applied to the 
model is equal to the mean quasi-static damage onset load determined from experiments. The 
load-displacement and the load-strain behavior computed from different FE models were 
compared to the corresponding experimental results. This global response was used to examine 
whether the FE models, the boundary conditions, the loads and the material properties used in 
the model yielded reasonable results. Strains were averaged from computed nodal point values 
over a length corresponding to the dimensions of the strain gages shown in Figure 2(a) [10]. 

The load versus displacement plot in Figure 12 shows that shell/3D model exhibits a 
stiffer behavior compared to results from the full three-dimensional model. The models also 
show a slightly stiffer response compared to the experiments (extensometer). This discrepancy 
may be explained by the fact that the material data used in the FE simulation originated from the 
literature. For a consistent simulation, material data should be taken from the batch of material 
that was used to manufacture the specimens. For comparison, results from 2D analysis were 
included in the plot of Figure 12 [22], Results from the plane-strain analysis indicated a stiffer 
behavior. This is caused by the constraints inherent to the plane-strain model, particularly in the 
±45° plies. The plane-stress model, which imposes the out-of-plane stresses to be zero and 
allows the strain to be the free parameter, exhibits a more compliant behavior. The influence of 
the assumptions made in developing two-dimensional finite element models on skin-stiffener 
debonding specimens was studied in detail in references [22], where geometrically nonlinear 
finite element analyses using two-dimensional plane-stress and plane-strain elements as well as 
three different generalized plane strain type approaches were performed. 

A comparison of measured strains at the surface of the flange (see Figure 2(a)) and 
computed results is plotted in Figure 13. In Figure 14, measured strains at the surface of the top 
45° skin ply near the flange tip (see Figure 2(a)) and computed surface strains are compared. For 
both locations, the strains calculated from the shell/3D model are in good agreement with the 
results from the full three-dimensional finite element model. The models show a slightly stiffer 
response compared to the strains measured during the experiments. As before, results from 2D 
analysis were included in the plots of Figure 13 and 14 for comparison [22], The plane-strain 
model showed a stiffer behavior yielding an upper bound while the plane-stress models were 
more compliant, yielding a lower bound. 

For the current shell/3D model, the desired reduction in computation time could not be 
achieved. This is partially caused by the fine three-dimensional mesh in the vicinity of the 
delaminated region, which was not replaced by shell elements. Additionally, the transition from 
the global shell element model to the local three-dimensional model was accomplished by multi- 
point constraints to enforce appropriate translations and rotations at the shell-solid interface, 
which increased the computation time. 
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4.1.2. Fractures Mechanics Analysis of the Specimen Subjected to Tension Loading 

During a series of nonlinear finite element analyses, strain energy release rates were 
computed at each front location for the loads applied in the experiments. The computed results 
were plotted versus the crack length in the center of the specimen (z=12.7 mm) as depicted in 
Figure 15(a) and across the width of the specimen for the longest delamination modeled (x=31.2 
mm) as depicted in Figure 15(b). 

The total energy release rates Gt=Gi+Gh+Gih, along the centerline of the specimen 
(z=12.7 mm) obtained from shell/3D analysis and three-dimensional simulations [16] are plotted 
in Figure 16 for comparison. After a small initial drop the computed total energy release rate 
increases sharply with delamination length, reaches a peak value and gradually decreases. The 
values from the shell/3D analyses are within 3.5% of the results obtained from the full three- 
dimensional simulations. Improved results may be obtained by enlarging the section modeled 
with solid elements. The values from plain strain and plane stress analysis are included in 
Figure 16 for comparison [11, 22], Qualitatively, all results follow the same trend. As before, the 
values from plane strain and plane stress analyses form upper and lower bounds except for very 
short delamination lengths. 

The variation of mixed mode ratio Gs /Gt with delamination length is shown in Figure 17. 
Here Gs denotes the sum of the in-plane shearing components Gn+Gm, and Gt denotes the total 
energy release rate Gi+Gn+Gm, where Gi is the opening mode. For two-dimensional analyses, 
where Gni= 0 , this definition is equal to the commonly used definition of the mixed mode ratio, 
Gh/Gt. For three-dimensional analysis, which also yields results for the scissoring mode Gni , the 
modified definition of Gs is introduced. The delamination initially starts with high shearing 
components, followed by a drop, which is equivalent to an increase in Gj. For longer 
delaminations a gradual increase in shearing components is observed. The values from the 
shell/3D analysis are in good agreement with the results obtained from the full three-dimensional 
simulations [16]. The values from plain strain and plane stress analysis were included in 
Figure 17 for comparison and follow the same trend [11]. The results from three-dimensional 
analysis in the center of the specimen width (z=12.7 mm) show a higher shearing component 
compared to the results from two-dimensional analysis. 

The distribution of the computed total energy release rate Gt across the width (z- 
coordinate) of the specimen is plotted in Figure 18 for the longest delamination modeled (x=31.2 
mm). Values at the free edge (z=0.0 mm and z=25.4 mm) have been excluded from the plot as 
the model is not fine enough to accurately calculate the energy release rates at the immediate 
edge. Across the width the computed total energy release rate gradually increases with z before it 
drops off near the free edge (z=25.4 mm). Results from reference [16] - obtained from a full 
three-dimensional model - have been included in the plot for comparison. The Gt values from 
the shell/3D analysis are within 6% of the results obtained from the full three-dimensional 
simulations. The distribution of the mixed-mode ratio Gs /Gt across the width (z-coordinate) of 
the specimen is plotted in Figure 19 for the same location (x=31.2 mm). Across the width the 
computed mixed mode ratio shows an almost constant value in the center of the specimen with 
increasing shearing components near the edges (z=0.0 mm and z=25.4 mm). The plotted mixed- 
mode ratios are in excellent agreement. 
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Results obtained from the shell/3D model with a refined mesh near the edges (see Figure 7) 
are shown in Figures 20 and 21. The distribution of the computed total energy release rate Gt 
across the width (z-coordinate) of the specimen is plotted in Figure 20 for the longest 
delamination modeled (x=3 1 .2 mm). Across the width the computed total energy release rate 
gradually increases with z before it first drops off near the free edges and then sharply increases 
in the zone of the mesh refinement in the immediate vicinity of the free edges (z=0.0 mm and 
z=25.4 mm). The model with the refined mesh near the edges yields additional data points where 
a steep increase in Gt is observed which could not be computed with the model used earlier. 
Computed energy release rates were omitted where the mesh transitions from coarse to fine, 
causing two gaps in the distribution at z ~ 2 mm and z ~ 23.5 mm. The data was omitted from the 
plot since the post processing algorithm used does not accurately calculate the energy release rate 
at nodal points which belong to elements with different width. The distribution of the mixed- 
mode ratio Gs /Gt across the width (z-coordinate) of the specimen is plotted in Figure 21 for the 
same location (x=31.2 mm). Across the width, the computed mixed mode ratio shows high 
shearing components near the edges, followed by a drop toward the center of the specimen, 
which is equivalent to an increase in opening mode I. The current model provided a smoother 
distribution across the width and yielded additional data points near the edges where a steep 
increase in mixed-mode ratio is observed, which could not be obtained from the coarse model 
discussed in the previous paragraph. The plotted mixed-mode ratios are in excellent agreement. 

A comparison of results obtained from the model with the refined mesh near the edges with 
results obtained from the model with the coarse mesh is provided Figure 22 for the full 3D 
analysis, and in Figure 23 for the shell/3D analysis. Results are in excellent agreement and 
complement each other when considering that the nodal point coordinates along the front are 
different for both models. 


4.2 Analysis of the Skin/Stringer Specimen Subjected to Three-Point Bending Load Using 
Shell/3D Modeling Technique 

The current analysis of the skin/stringer specimen subjected to three-point bending load 
supplements the study in reference [16]. The focus of the current analysis is the comparison of 
mixed-mode energy release rates from reference [16] with those obtained using a new shell to 
solid coupling option offered by ABAQUS®6.3 [18]. 


4.2.1. Global Response of the Specimen Subjected to Three-Point Bending Load 

The global response of the specimens was computed at the mean quasi-static damage 
onset load determined from experiments. The load-displacement and the load-strain behavior 
computed from different FE models were compared to the corresponding experimental results. 
This global response was used to examine whether the FE models, the boundary conditions, the 
loads and the material properties used in the models yielded reasonable results. Strains were 
averaged from computed nodal point values over a length corresponding to the dimensions of the 
strain gages shown in Figure 2(a) [10]. 
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The load versus displacement data is plotted in Figure 24. Results from the current 
shell/3D analysis, which used the new surface based option in ABAQUS® to connect the global 
shell model to the local solid model, are labeled (*SHELL TO SOLID). For comparison, results 
from reference [16] where the shell and solid regions were connected using multi-point 
constraints (*MPC) in ABAQUS® were included in the figure. Both shell/3D models are in 
excellent agreement and exhibit a stiffer behavior compared to results from the full three- 
dimensional model. The models also predict a slightly stiffer response compared to the 
experiments (DCDT). This discrepancy may be explained by the fact that the material data used 
in the FE simulation originated from the literature. For a consistent simulation, material data 
should be taken from the batch of material that was used to manufacture the specimens. For 
comparison, results from 2D analysis were included in the plot of Figure 24 [22], Results from 
the plane-strain analysis indicated a stiffer behavior. This is caused by the constraints inherent to 
the plane-strain model, particularly in the ±45° plies. The plane-stress model which imposes the 
out of plane stresses to be zero and allows the strain to be the free parameter exhibits a more 
compliant behavior compared to results from full three-dimensional analysis. 

A comparison of measured strains at the surface of the flange (see Figure 2(a)) and 
computed results is provided in Figure 25. In Figure 26, measured strains at the surface of the top 
45° skin ply near the flange tip (see Figure 2(a)) and computed surface strains were compared. 
For both locations, the strains calculated from the both shell/3D models are in excellent 
agreement with the results from the full three-dimensional finite element model. The models 
show a slightly stiffer response compared to the experiments. As before, results from 2D 
analysis were included in the plots of Figure 25 and 26 for comparison [22], The plane-strain 
model predicted a stiffer behavior, yielding an upper bound, while the plane-stress models were 
more compliant, yielding a lower bound. 


4.2.2. Fractures Mechanics Analysis of the Specimen Subjected to Three-Point Bending 
Load 

The distribution of the computed total energy release rate Gt across the width (z- 
coordinate) of the specimen is plotted in Figure 27 for the longest delamination modeled (x=31.2 
mm) as shown in Figure 15(b). Across the width the computed total energy release rate gradually 
increases with z before it drops off near the free edge (z = 25.4 mm). Results from references 
[22] - obtained from the full three-dimensional models - have been included in the plot for 
comparison. The Gt values from shell/3D analyses are in excellent agreement and are about 6% 
lower than the results obtained from the full three-dimensional simulations. The distribution of 
the mixed-mode ratio Gs /Gt across the width (z-coordinate) of the specimen is plotted in 
Figure 28 for the same location (x=3 1.2 mm). Across the width the computed mixed mode ratios 
are in excellent agreement and show an almost constant value with increasing shearing 
components near the edges (z=0.0 mm andz=25.4 mm). 

For the longest delamination modeled (x=3 1 .2 mm), the distribution of the computed 
total energy release rate Gt across the width (z-coordinate) of the specimen is plotted in 
Figure 29. Results were obtained from the shell/3D model with a refined mesh near the edges as 
shown in Figure 7. Across the width, the computed total energy release rate gradually increases 
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before it first drops off near the free edges and then sharply increases in the zone of the mesh 
refinement in the immediate vicinity of the free edges (z=0.0 mm and z=25.4 mm). The model 
with the refined mesh near the edges yielded additional data points where a steep increase in Gt 
is observed which could not be computed with the less refined model (Figure 27). The shell/3D 
yields results are in good agreement with energy release rates obtained from a full 3D model. 
Computed energy release rates were omitted where the mesh transitions from coarse to fine 
causing two gaps in the distribution at z ~ 2 mm and z ~ 23.5 mm. The data was omitted from the 
plot since the post processing algorithm used does not accurately calculate the energy release rate 
at nodal points which belong to elements with different width. The distribution of the mixed- 
mode ratio Gs /Gt across the width (z-coordinate) of the specimen is plotted in Figure 30 for the 
same location (x=31.2 mm). Across the width the computed mixed mode ratio is nearly constant 
and dominated by opening mode I, but shows high shearing components near the edges, which is 
caused by the overclosure of the delaminated surface near the edges where the crack opening 
mode I disappears. The current model provided a smoother distribution across the width and 
yielded additional data points near the edges where a steep increase in mixed-mode ratio Gs /Gt 
is observed which could not be computed with the less refined model (Figure 28). 


4.3. Analysis of the Skin/Stringer Specimen Subjected to Tension Load Using Submodeling 
Technique 

4.3.1. Computed Displacements and Strains for the Specimen Subjected to Tension 
Loading 

Initially, the global response of the specimens was computed at the mean quasi-static 
damage onset load determined from experiments. The load-displacement and the load-strain 
behavior computed from different global models were compared to the corresponding 
experimental results. This global response was used to examine whether the global models, the 
boundary conditions, the loads and the material properties used in the model produced 
reasonable results. Strains were averaged from computed nodal point values over a length 
corresponding to the dimensions of the strain gages shown in Figure 2(a). 

The load versus displacement plot is shown in Figure 31. Results show that both global 
shell models are in excellent agreement and exhibit a slightly stiffer behavior compared to results 
from the full three-dimensional global model and results from full 3D analysis published in 
reference [16]. All FE models show a stiffer response compared to the experiments 
(extensometer). This discrepancy may be explained by the fact that the material data used in the 
FE simulation originated from the literature. For a consistent simulation, material data should be 
taken from the batch of material that was used to manufacture the specimens. Due to matching 
results, the global model made of S4 elements was not studied further. Further investigations 
focussed on the global model made of S8R elements. 

A comparison of measured strains at the surface of the flange (see Figure 2(a)) and 
computed results is shown in Figure 32. In Figure 33, measured strains at the surface of the top 
45° skin ply near the flange tip (see Figure 2(a)) and computed surface strains were compared. 
For both locations, the strains calculated from the global shell model are in excellent agreement 
with the results from the full three-dimensional finite element model. The models show a slightly 
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stiffer response compared to the strains measured during the experiments. Based on the 
matching displacements discussed above, the strains from the global solid model were assumed 
to be identical to the strains obtained from the full 3D analysis of reference [16] and hence were 
not included in the figures. 

The global models made of shell and solid elements appeared to model the displacement 
and strain behavior of the skin/stringer specimen subjected to tension with sufficient accuracy. 
Therefore the results from both global models were used as boundary conditions for the local 
submodels. 


4.3.2. Computed Mixed-Mode Strain Energy Release Rates for the Specimen Subjected to 
Tension Loading 

During a series of nonlinear finite element analyses, strain energy release rates were 
computed at each front location for the loads applied in the experiments. In subsequent analyses 
results from the global shell model as shown in Figure 9(a) or results from the global solid model 
as shown in Figure 11(b) were used as boundary conditions for the submodels. For the global 
shell model, the computed strain energy release rate distributions across the width of the 
specimen (atx=31 .2 mm, as shown in Figure 15(b)) are plotted in Figure 34 for different sizes of 
the submodels as shown in Figures 9 and 10. For the global solid model the computed results are 
depicted in Figure 35. Values at the free edge (z= 0.0 mm and z=25.4 mm) have been excluded 
from the plots as the model was not fine enough to accurately model the influence of the free 
edges on the distribution of the energy release rates. With increasing length of the local three- 
dimensional submodel, an increase in computed strain energy release rate is observed. However, 
even for a submodel which extended six specimen thicknesses in front and behind the damaged 
region, the results did not converge to the values obtained from reference [16] which were 
included in Figures 34 and 35 for comparison. For the global shell model, the distribution of the 
mixed-mode ratio Gs /Gt across the width (z-coordinate) of the specimen is plotted in Figure 36 
for the same location (x=3 1 .2 mm). For the global solid model the computed results are depicted 
in Figure 37. With increasing length of the local three-dimensional submodel, computed results 
converge to mixed-mode ratios obtained from reference [16] which were included in Figures 36 
and 37 for comparison. 


4.4. Analysis of the Skin/Stringer Specimen Subjected to Three-Point Bending Load Using 
Submodeling 

4.4.1. Computed Displacements and Strains for the Specimen Subjected to Three-Point 
Bending Load 

The global response of the specimen was computed at the mean quasi-static damage 
onset load determined from experiments. The load-displacement and the load-strain behavior 
computed from different FE models were compared to the corresponding experimental results. 
This global response was used to examine whether the FE models, the boundary conditions, the 
loads and the material properties used in the models yielded reasonable results. Strains were 
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averaged from computed nodal point values over a length corresponding to the dimensions of the 
strain gages shown in Figure 2(a). 

The load versus displacement plot is shown in Figure 38. Results show that both global 
shell models are in excellent agreement and exhibit a stiffer behavior compared to results from 
the full three-dimensional global model and the full 3D analysis from reference [16]. 

A comparison of measured strains at the surface of the flange (see Figure 2(a)) and 
computed results is shown in Figure 39. In Figure 40, measured strains at the surface of the top 
45° skin ply near the flange tip (see Figure 2(a)) and computed surface strains were compared. 
For both locations, the strains calculated from both global shell models are in excellent 
agreement with the results from the full three-dimensional finite element model from reference 
[16]. The models show a slightly stiffer response compared to the strains measured during the 
experiments. Based of the matching displacements discussed above, the strains from the global 
solid model were assumed to be identical to the strains obtained from the full 3D analysis from 
reference [16] and hence were not included in the figures. Due to matching results the global 
model made of S4 elements was not studied further and the further investigations focused on the 
global model made of S8R elements. 

The global models made of shell and solid elements appeared to model the displacement 
and strain behavior of the skin/stringer specimen subjected to three-point bending sufficiently 
accurate. Therefore the results from both global models were used as boundary conditions for the 
local submodels. 


4.3.2. Computed Mixed-Mode Strain Energy Release Rates for the Specimen Subjected to 
Three-Point Bending Load 

During a series of nonlinear finite element analyses, strain energy release rates were 
computed at each front location for the loads applied in the three-point bending test. In 
subsequent analyses results from the global shell model as shown in Figure 9(a) or results from a 
global solid model as shown in Figure 11(b) were used as boundary conditions of the submodels. 
For the global shell model, the computed strain energy release rate distributions across the width 
of the specimen (at x=31.2 mm, as shown in Figure 15(b)) are plotted in Figure 41 for different 
sizes of the submodels as shown in Figures 9 and 10. For the global solid model, computed total 
strain energy release rate distributions across the width (z-coordinate) of the specimen are 
depicted in Figure 42 at the same location (x=31.2 mm) for different sizes of the submodels. 
Values at the free edge (z=0.0 mm and z=25.4 mm) have been excluded from the plots, as the 
model was not fine enough to accurately model the influence of the free edges on the distribution 
of the energy release rates. With increasing length of the local three-dimensional submodel, an 
increase of the computed total strain energy release rate is observed. However, even for a 
submodel which extends six specimen thickness in front and behind the damaged region the 
results did not converge to the solution obtained from reference [16], which was included in 
Figures 41 and 42 for comparison. For the local shell model, the distribution of the mixed-mode 
ratio Gs /Gt across the width (z-coordinate) of the specimen is plotted in Figure 43 for the same 
location (x=31.2 mm). For the global solid model the computed results are depicted in Figure 44. 
With increasing length of the local three-dimensional submodel, computed results converge to 
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the solution obtained from reference [16], which was included in Figures 43 and 44 for 
comparison. 


5. SUMMARY AND CONCLUDING REMARKS 

The application of the shell/3D modeling technique for the simulation of skin/stringer 
debond specimen subjected to tension and three-point bending was demonstrated. The study 
extended the application of this technique beyond the simulation of simple specimens - such as 
Double Cantilever Beam (DCB), End Notched Flexure (ENF) and Single Leg Bending (SLB) 
specimens - where the delamination is located between unidirectional plies in the mid-plane of 
the specimen to a coupon type specimen where is the delamination is located between plies of 
different orientation. The global structure was modeled with shell elements. A three-dimensional 
model, extending to about three specimen thicknesses on either side of the delamination front 
was used locally. Delaminations of various lengths were discretely modeled at the locations 
where delaminations were observed during previous experiments. 

Computed strains at the center of the stringer flange were compared with results obtained 
from full three-dimensional finite element analysis and experimental data. Results from analyses 
and experiment were in good agreement. Calculated total strain energy release rates and mixed 
mode ratios obtained form shell/3D simulations also were in good agreement with results 
obtained from full solid models. Both options offered by ABAQUS® to connect the shell and 
solid regions yielded almost identical results. The concurrence of the results demonstrated the 
effectiveness of the shell/3D modeling technique for the investigation of delamination onset from 
an initial crack. For comparison, results from plane stress and plane strain models provided 
displacements, flange strains, as well as energy release rates, which formed an upper and lower 
bound of the results obtained from full three-dimensional or shell/3D simulations. 

The current study provided an additional verification step for the shell/3D modeling 
technique prior to its application to large, full-scale stringer stiffened panels based on the 
assumption that the failure observed in the skin/stringer specimen is identical to the failure 
observed in the full-scale panels where the delamination is located between plies of different 
orientation. For the current shell/3D model, however, the desired reduction in computation time 
could not be achieved. This is partially caused by the fine three-dimensional mesh in the vicinity 
of the delaminated region, which was not replaced by shell elements. Additionally, the transition 
from the global shell element model to the local three-dimensional model was accomplished by 
multi-point constraints to enforce appropriate translations and rotations at the shell-solid 
interface, which increased the computation time. As the model size increases, however, gains in 
efficiency will be more significant. For large, built-up composite structures modeled with plate 
elements, the shell/3D modeling technique may be used for reducing the model size, compared to 
a full three-dimensional model, since only a relatively small section in the vicinity of the 
delamination front needs to be modeled with solid elements. Efficiency gains are expected to be 
larger for larger models. 

The application of the submodeling technique for the simulation of skin/stringer debond 
in a specimen subjected to tension and three-point bending was also studied. Global models 
made of shell elements and solid elements were used. Solid elements were used for local 
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submodels, which extended between three and six specimen thicknesses on either side of the 
delamination front to model the details of the damaged section. Computed total strain energy 
release rates did not converge to results obtained from full solid models. Computed mixed-mode 
ratios however, appeared to slowly converge to the results obtained from full solid models. 
Based on the current results presented, the submodeling technique cannot be recommended for 
large, built-up composite structures. 
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Table 1. 


Material Properties 


IM7/8552 Unidirectional Graphite/Epoxy Prepreg 

E\i = 161.0 GPa 

^22 = 11.38 GPa 

£33 = 11.38 GPa 

vi 2 = 0.32 

vi 3 = 0.32 

V23 = 0.45 

Gi 2 = 5.17 GPa 

Gn = 5.17 GPa 

G 23 = 3.92 GPa 

IM7/8552 Graphite/Epoxy Plain Weave Fabric 

En = 71.7 GPa 

£22 = 71.7 GPa 

£33 = 10.3 GPa 

vi 2 = 0.04 

vi 3 = 0.35 

V23 = 0.35 

Gi 2 = 4.48 GPa 

Gi 3 = 4.14 GPa 

G 23 = 4.14 GPa 

Grade 5 FM300 Adhesive 

£=1.72 GPa 

v= 0.3 

(assumed isotropic) 


The material properties are given with reference to the ply coordinate axes where index 1 1 denotes 
the ply principal axis that coincides with the direction of maximum in-plane Young’s modulus 
(fiber direction). Index 22 denotes the direction transverse to the fiber in the plane of the lamina 
and index 33 the direction perpendicular to the plane of the lamina. 
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u, v: displacements 


y.v 


Lx,U 



u=v=0 I Q=427.8 N v=0 

(a) Deformed two-dimensional model of skin/ 'stringer specimen (101.6 mm x 25.4 mm) [10]. 

v=0 for 0< z <25.4 mm 
w=0 at z=12.7 mm 



v=0 for 0< z <25.4 mm 
u=w=0 at z=12.7 mm 


Q=427.8 N 


z,w 


undeformed outline 
u, v, w: displacements 


(b) Deformed three-dimensional model of skin/ 'stringer specimen (101.6 mm x 25.4 mm) [11 ]. 

global shell element model of 

unfractured section 

local 3D model of 

delaminated section 

v=0 for 0< z <25.4 mm 
w=0 at z=12.7 mm 



Q=427.8 N 


undeformed outline 
u, v, w: displacements 


v=0 for 0< z <25.4 mm 

z,w u=w=0 at z=12.7 mm 

(c) Deformed shell/3D model of skin/ stringer specimen (101.6 mm x 25.4 mm) [ 16]. 



Figure 1. Finite Element Modeling. 
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(a). Specimen configuration. 
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Q= 427.6 N 


101 .6 mm 

(c). Three-point bending case. 

Figure 2. Specimen configuration and test set-up. 
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(a). Typical damage patterns observed at specimen corners 1 and 4. 
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(b) . Typical damage patterns observed at specimen corners 2 and 3. 
Figure 3. Typical damage patterns observed in skin/ stringer specimen [10]. 
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Shell to solid 
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(a). Outline ofshell/3D model. 


u, v, w: displacements 


global shell element model of 
unfractured section (S8R elements) 



















local 3D model of 
delaminated section 
(C3D20R elements) 

u, v, w: displacements 


(b). Shell/3D finite element model. 



(c). Detail of undeformed local 3D model and shell to 3D interface region. 
Figure 5. Shell/3D model of skin/flange specimen. 
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(a). Deformed shell/3 D model of tension specimen with load and boundary conditions. 


y.v 

x,u 
z,w 

(b). Deformed shell/3D model of bending specimen with load and boundary conditions. 





(c). Detail of deformed model of damaged specimen. 
Figure 6. Deformed shell/3 D model of skin/flange specimen. 
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(a). Deformed shell/3D model with load and boundary conditions. 



(b). Detail of undeformed local 3D model and shell to 3D interface region. 



(c). Detail of deformed model of damaged specimen. 

Figure 7. Shell/3 D model of skin/flange specimen with fine mesh near the edge. 
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entire specimen modeled with shell elements 


y.v 


t^x,u 


z,w 



physical location of submodel 


u, v, w: displacements 


(a). Outline of submodelling technique using shell elements for global model 
and solid elements for the submodel. 




Figure 8. Submodelling technique using a global model made of shell elemets. 
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z,w 


v=0 for 0< z <25.4 mm 
w=0 at z=12.7 mm 


Q=427.8 N 


v=0 for 0< z <25.4 mm 
u=w=0 at z=12.7 mm 


— ■ - undeformed outline 
u, v, w: displacements 


(a). Deformed global shell model of bending specimen with load and boundary conditions. 



(b). Deformed global shell model of tension specimen with load and boundary conditions 



s: skin thickness 
f: flange thickness 

maximum length of delamination 
modeled 

matrix crack and origin of delamination 
(c). Submodel made of solids extending about three specimen thicknesses beyond damage. 
Figure 9. Finite element models used for submodeling technique. 
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y.v 

k 


x,u 



z,w 

u, v, w: displacements 


s: skin thickness 
f: flange thickness 


maximum length of delamination modeled 
matrix crack and origin of delamination 


(a). Submodel made of solids extending about four specimen thicknesses beyond damage. 



s: skin thickness 
f: flange thickness 


maximum length of delamination modeled 
matrix crack and origin of delamination 


u, v, w: displacements 

(b). Submodel made of solids extending about five specimen thicknesses beyond damage. 



z,w 

u, v, w: displacements 


s: skin thickness 
f: flange thickness 


maximum length of delamination modeled 
matrix crack and origin of delamination 


(c). Submodel made of solid elements extending about six specimen thicknesses beyond damage. 
Figure 10. Submodels made of solid elements. 
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solid elements used for global model 


y.v 
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x,u 


z,w 



physical location of submodel 


u, v, w: displacements 

(a). Outline of submodelling technique using solid elements for global model and submodel. 



(b). Global model of skin/flange specimen modeled with solid elements. 



(c). Outline of submodel with discretely modeled damage. 

Figure 11. Submodelling technique using a global model made of solid elemets. 
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Figure 13. Flange strain-load plots for tension tests. 



Load, kN 

Figure 14. Typical skin strain-load plot for tension tests. 
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(a). XY-plane at z=12.5 mm for display of results along the length in the center of the specimen. 



(b). YZ-plane at x=3 1.1 mm for display of results across the width of the specimen. 
Figure 15. Section planes through the specimen 
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Figure 16. Computed total energy release rate for delamination growing from matrix crack 
in skin top 45°/-45°ply interface for tension test. 
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Figure 17. Computed mixed mode ratio G/G^for delamination growing 
from matrix crack in skin top 45°/-45°ply interface for tension test. 


31 



coordinate z, mm 

Figure 18. Computed total energy release rate across the width of the specimen 
at x= 31.2 mm for tension test (coarse mesh). 
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Figure 19. Computed mixed mode ratio across the width of the specimen 
at x= 31.2 mm for tension test ( coarse mesh). 
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coordinate z, mm 

Figure 20. Computed total energy release rate across the width of the specimen 
at x= 31.2 mm for tension test ( refined mesh). 



coordinate z, mm 

Figure 21. Computed mixed-mode ratio across the width of the specimen 
at x= 31.2 mm for tension test (refined mesh). 
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coordinate z, mm 

Figure 22. Computed total energy release rate across the width of the specimen 
at x= 31.2 mm for tension test (full 3D analysis). 



coordinate z, mm 

Figure 23. Computed total energy release rate across the width of the specimen 
at x= 31.2 mm for tension test (shell/3D analysis ). 
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Figure 24. Load-displacement plots for three-point bending tests. 


35 


4000 


3000 


Strain 
x 10' 6 


2000 




Flange 

debonding 


»« 


1000 


2D plain strain elements 
2D plain stress elements 
" • full 3D analysis 

shell/3D analysis (*MPC) 

shell/3D analysis (*SHELL TO SOLID) 


0*_* — * — * — * — 1 — * — * — * — * — 1 — * — * — * — * — 1 — * — * — * — * — 1 — * — * — * — * — l- 

0 100 200 300 400 500 


Load, N 


Figure 25. Flange strain-load plots for three-point bending tests. 
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Figure 27. Computed total energy release rate across the width of the specimen 
at x= 31.2 mm for three-point bending test ( coarse mesh). 
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Figure 28. Computed mixed mode ratio across the width of the specimen 
at x= 31.2 mm for three-point bending test ( coarse mesh). 
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Figure 29. Computed total energy release rate across the width of the specimen 
at x= 31.2 mm for three-point bending test (refined mesh). 
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Figure 30. Computed mixed mode ratio across the width of the specimen 
at x= 31.2 mm for three-point bending test (refined mesh). 
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Figure 32. Flange strain-load plots for tension tests. 
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Figure 34. Computed total energy release rate across the width of the specimen 

at x= 31.2 mm for tension test. 
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Figure 35. Computed total energy release rate across the width of the specimen 

at x= 31.2 mm for tension test. 
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Figure 36. Computed mixed mode ratio across the width of the specimen 
at x= 31.2 mm for tension test. 
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Figure 37. Computed mixed mode ratio across the width of the specimen 
at x= 31.2 mm for tension test. 
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Figure 38. Load-displacement plots for three-point bending tests. 
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Figure 39. Flange strain-load plots for three-point bending tests. 
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Figure 41. Computed total energy release rate across the width of the specimen 
at x= 31.2 mm for three-point bending test. 
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Figure 42. Computed total energy release rate across the width of the specimen 
at x= 31.2 mm for three-point bending test. 
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Figure 43. Computed mixed mode ratio across the width of the specimen 
at x= 31.2 mm for three-point bending test. 
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Figure 44. Computed mixed mode ratio across the width of the specimen 
at x= 31.2 mm for three-point bending test. 
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